from sympy import symbols, sqrt
from sympy.simplify import ratsimp
from sympy.simplify.cse_main import cse

vx1, vy1, vz1 = symbols('vx1 vy1 vz1', real=True)
vx2, vy2, vz2 = symbols('vx2 vy2 vz2', real=True)
vx3, vy3, vz3 = symbols('vx3 vy3 vz3', real=True)
vx4, vy4, vz4 = symbols('vx4 vy4 vz4', real=True)

xa = vx2-vx3
ya = vy2-vy3
za = vz2-vz3
x1 = vx1-vx2
y1 = vy1-vy2
z1 = vz1-vz2
x2 = vx4-vx3
y2 = vy4-vy3
z2 = vz4-vz3

nx1 = y1 * za - z1 * ya
ny1 = z1 * xa - x1 * za
nz1 = x1 * ya - y1 * xa

nx2 = y2 * za - z2 * ya
ny2 = z2 * xa - x2 * za
nz2 = x2 * ya - y2 * xa

vn1xvaX = ny1*za - nz1*ya
vn1xvaY = nz1*xa - nx1*za
vn1xvaZ = nx1*ya - ny1*xa

vn1xv1X = ny1*z1 - nz1*y1
vn1xv1Y = nz1*x1 - nx1*z1
vn1xv1Z = nx1*y1 - ny1*x1

vn1xv2X = ny1*z2 - nz1*y2
vn1xv2Y = nz1*x2 - nx1*z2
vn1xv2Z = nx1*y2 - ny1*x2

vn2xvaX = ny2*za - nz2*ya
vn2xvaY = nz2*xa - nx2*za
vn2xvaZ = nx2*ya - ny2*xa

vn2xv2X = ny2*z2 - nz2*y2
vn2xv2Y = nz2*x2 - nx2*z2
vn2xv2Z = nx2*y2 - ny2*x2

vn2xv1X = ny2*z1 - nz2*y1
vn2xv1Y = nz2*x1 - nx2*z1
vn2xv1Z = nx2*y1 - ny2*x1


r1 = sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
r2 = sqrt(nx2*nx2+ny2*ny2+nz2*nz2)
r12 = sqrt((nx1*nx1 + ny1*ny1 + nz1*nz1)*(nx2*nx2 + ny2*ny2 + nz2*nz2))

dotN12 = (nx1*nx2+ny1*ny2+nz1*nz2)
cost = dotN12/sqrt((nx1*nx1+ny1*ny1+nz1*nz1)*(nx2*nx2+ny2*ny2+nz2*nz2))

symbols = {r12:'r12', dotN12:'dotN12', r1:'r1', r2:'r2', 
    vn1xvaX:'vn1xvaX', vn1xvaY:'vn1xvaY', vn1xvaZ:'vn1xvaZ',
    vn1xv1X:'vn1xv1X', vn1xv1Y:'vn1xv1Y', vn1xv1Z:'vn1xv1Z',
    vn2xvaX:'vn2xvaX', vn2xvaY:'vn2xvaY', vn2xvaZ:'vn2xvaZ',
    vn2xv2X:'vn2xv2X', vn2xv2Y:'vn2xv2Y', vn2xv2Z:'vn2xv2Z',
    (nz2*ya - ny2*za):'-vn2xvaX', (nx2*za - nz2*xa):'-vn2xvaY', (ny2*xa - nx2*ya):'-vn2xvaZ',
    (nz1*ya - ny1*za):'-vn1xvaX', (nx1*za - nz1*xa):'-vn1xvaY', (ny1*xa - nx1*ya):'-vn1xvaZ',
    (nz1*y1 - ny1*z1):'-vn1xv1X', (nx1*z1 - nz1*x1):'-vn1xv1Y', (ny1*x1 - nx1*y1):'-vn1xv1Z',
    (nz2*y2 - ny2*z2):'-vn2xv2X', (nx2*z2 - nz2*x2):'-vn2xv2Y', (ny2*x2 - nx2*y2):'-vn2xv2Z',
    (ny1*z2 + ny2*z1 - nz1*y2 - nz2*y1):'(vn1xv2X + vn2xv1X)',
    (nz1*x2 + nz2*x1 - nx1*z2 - nx2*z1):'(vn1xv2Y + vn2xv1Y)',
    (nx1*y2 + nx2*y1 - ny1*x2 - ny2*x1):'(vn1xv2Z + vn2xv1Z)',
    nx1:'nx1', ny1:'ny1', nz1:'nz1', 
    nx2:'nx2', nz2:'nz2', nz2:'nz2',
    xa*z2 - x2*za:'ny2',
    (-3*nz1*ya + 3*ny1*za):'3*vn1xvaX',  (-3*nx1*za + 3*nz1*xa):'3*vn1xvaY',  (-3*ny1*xa + 3*nx1*ya):'3*vn1xvaZ',
    (-3*ny1*z1 + 3*nz1*y1):'-3*vn1xv1X', (-3*nz1*x1 + 3*nx1*z1):'-3*vn1xv1Y', (-3*nx1*y1 + 3*ny1*x1):'-3*vn1xv1Z', 
    (-3*ny2*z2 + 3*nz2*y2):'-3*vn2xv2X', (-3*nz2*x2 + 3*nx2*z2):'-3*vn2xv2Y', (-3*nx2*y2 + 3*ny2*x2):'-3*vn2xv2Z', 
}

old = dotN12*(-ya**2 - za**2)/(r1**3*r2) + dotN12*vn1xvaX*3*vn1xvaX/(r1**5*r2) - 2*vn1xvaX*vn2xvaX/(r1**3*r2)

print "dx2[0][0][0] =",        cost.diff(vx1, vx1).subs(symbols)
print "dx2[0][0][1] =",        cost.diff(vx1, vy1).subs(symbols)
print (cost.diff(vx1, vy1) - (dotN12*(xa*ya +          3*vn1xvaX*vn1xvaY/(r1**2)) - vn1xvaX*vn2xvaY - vn1xvaY*vn2xvaX )/(r1**3*r2)).expand()
"""
print "dx2[0][0][2] =",        cost.diff(vx1, vz1).subs(symbols)
print "dx2[0][1][0] =",        cost.diff(vy1, vx1).subs(symbols)
print "dx2[0][1][1] =",        cost.diff(vy1, vy1).subs(symbols)
print "dx2[0][1][2] =",        cost.diff(vy1, vz1).subs(symbols)
print "dx2[0][2][0] =",        cost.diff(vz1, vx1).subs(symbols)
print "dx2[0][2][1] =",        cost.diff(vz1, vy1).subs(symbols)
print "dx2[0][2][2] =",        cost.diff(vz1, vz1).subs(symbols)
        
print "dx2[1][0][0] =",        cost.diff(vx1, vx2).subs(symbols)
print "dx2[1][0][1] =",        cost.diff(vx1, vy2).subs(symbols)
print "dx2[1][0][2] =",        cost.diff(vx1, vz2).subs(symbols)
print "dx2[1][1][0] =",        cost.diff(vy1, vx2).subs(symbols)
print "dx2[1][1][1] =",        cost.diff(vy1, vy2).subs(symbols)
print "dx2[1][1][2] =",        cost.diff(vy1, vz2).subs(symbols)
print "dx2[1][2][0] =",        cost.diff(vz1, vx2).subs(symbols)
print "dx2[1][2][1] =",        cost.diff(vz1, vy2).subs(symbols)
print "dx2[1][2][2] =",        cost.diff(vz1, vz2).subs(symbols)
     
print "dx2[2][0][0] =",        cost.diff(vx1, vx3).subs(symbols)
print "dx2[2][0][1] =",        cost.diff(vx1, vy3).subs(symbols)
print "dx2[2][0][2] =",        cost.diff(vx1, vz3).subs(symbols)
print "dx2[2][1][0] =",        cost.diff(vy1, vx3).subs(symbols)
print "dx2[2][1][1] =",        cost.diff(vy1, vy3).subs(symbols)
print "dx2[2][1][2] =",        cost.diff(vy1, vz3).subs(symbols)
print "dx2[2][2][0] =",        cost.diff(vz1, vx3).subs(symbols)
print "dx2[2][2][1] =",        cost.diff(vz1, vy3).subs(symbols)
print "dx2[2][2][2] =",        cost.diff(vz1, vz3).subs(symbols)
        
print "dx2[3][0][0] =",        cost.diff(vx1, vx4).subs(symbols)
print "dx2[3][0][1] =",        cost.diff(vx1, vy4).subs(symbols)
print "dx2[3][0][2] =",        cost.diff(vx1, vz4).subs(symbols)
print "dx2[3][1][0] =",        cost.diff(vy1, vx4).subs(symbols)
print "dx2[3][1][1] =",        cost.diff(vy1, vy4).subs(symbols)
print "dx2[3][1][2] =",        cost.diff(vy1, vz4).subs(symbols)
print "dx2[3][2][0] =",        cost.diff(vz1, vx4).subs(symbols)
print "dx2[3][2][1] =",        cost.diff(vz1, vy4).subs(symbols)
print "dx2[3][2][2] =",        cost.diff(vz1, vz4).subs(symbols)
     
print "dx2[4][0][0] =",        cost.diff(vx2, vx1).subs(symbols)
print "dx2[4][0][1] =",        cost.diff(vx2, vy1).subs(symbols)
print "dx2[4][0][2] =",        cost.diff(vx2, vz1).subs(symbols)
print "dx2[4][1][0] =",        cost.diff(vy2, vx1).subs(symbols)
print "dx2[4][1][1] =",        cost.diff(vy2, vy1).subs(symbols)
print "dx2[4][1][2] =",        cost.diff(vy2, vz1).subs(symbols)
print "dx2[4][2][0] =",        cost.diff(vz2, vx1).subs(symbols)
print "dx2[4][2][1] =",        cost.diff(vz2, vy1).subs(symbols)
print "dx2[4][2][2] =",        cost.diff(vz2, vz1).subs(symbols)
        
print "dx2[5][0][0] =",        cost.diff(vx2, vx2).subs(symbols)
print "dx2[5][0][1] =",        cost.diff(vx2, vy2).subs(symbols)
print "dx2[5][0][2] =",        cost.diff(vx2, vz2).subs(symbols)
print "dx2[5][1][0] =",        cost.diff(vy2, vx2).subs(symbols)
print "dx2[5][1][1] =",        cost.diff(vy2, vy2).subs(symbols)
print "dx2[5][1][2] =",        cost.diff(vy2, vz2).subs(symbols)
print "dx2[5][2][0] =",        cost.diff(vz2, vx2).subs(symbols)
print "dx2[5][2][1] =",        cost.diff(vz2, vy2).subs(symbols)
print "dx2[5][2][2] =",        cost.diff(vz2, vz2).subs(symbols)
        
print "dx2[6][0][0] =",        cost.diff(vx2, vx3).subs(symbols)
print "dx2[6][0][1] =",        cost.diff(vx2, vy3).subs(symbols)
print "dx2[6][0][2] =",        cost.diff(vx2, vz3).subs(symbols)
print "dx2[6][1][0] =",        cost.diff(vy2, vx3).subs(symbols)
print "dx2[6][1][1] =",        cost.diff(vy2, vy3).subs(symbols)
print "dx2[6][1][2] =",        cost.diff(vy2, vz3).subs(symbols)
print "dx2[6][2][0] =",        cost.diff(vz2, vx3).subs(symbols)
print "dx2[6][2][1] =",        cost.diff(vz2, vy3).subs(symbols)
print "dx2[6][2][2] =",        cost.diff(vz2, vz3).subs(symbols)
        
print "dx2[7][0][0] =",        cost.diff(vx2, vx4).subs(symbols)
print "dx2[7][0][1] =",        cost.diff(vx2, vy4).subs(symbols)
print "dx2[7][0][2] =",        cost.diff(vx2, vz4).subs(symbols)
print "dx2[7][1][0] =",        cost.diff(vy2, vx4).subs(symbols)
print "dx2[7][1][1] =",        cost.diff(vy2, vy4).subs(symbols)
print "dx2[7][1][2] =",        cost.diff(vy2, vz4).subs(symbols)
print "dx2[7][2][0] =",        cost.diff(vz2, vx4).subs(symbols)
print "dx2[7][2][1] =",        cost.diff(vz2, vy4).subs(symbols)
print "dx2[7][2][2] =",        cost.diff(vz2, vz4).subs(symbols)
        
print "dx2[8][0][0] =",        cost.diff(vx3, vx1).subs(symbols)
print "dx2[8][0][1] =",        cost.diff(vx3, vy1).subs(symbols)
print "dx2[8][0][2] =",        cost.diff(vx3, vz1).subs(symbols)
print "dx2[8][1][0] =",        cost.diff(vy3, vx1).subs(symbols)
print "dx2[8][1][1] =",        cost.diff(vy3, vy1).subs(symbols)
print "dx2[8][1][2] =",        cost.diff(vy3, vz1).subs(symbols)
print "dx2[8][2][0] =",        cost.diff(vz3, vx1).subs(symbols)
print "dx2[8][2][1] =",        cost.diff(vz3, vy1).subs(symbols)
print "dx2[8][2][2] =",        cost.diff(vz3, vz1).subs(symbols)
        
print "dx2[9][0][0] =",        cost.diff(vx3, vx2).subs(symbols)
print "dx2[9][0][1] =",        cost.diff(vx3, vy2).subs(symbols)
print "dx2[9][0][2] =",        cost.diff(vx3, vz2).subs(symbols)
print "dx2[9][1][0] =",        cost.diff(vy3, vx2).subs(symbols)
print "dx2[9][1][1] =",        cost.diff(vy3, vy2).subs(symbols)
print "dx2[9][1][2] =",        cost.diff(vy3, vz2).subs(symbols)
print "dx2[9][2][0] =",        cost.diff(vz3, vx2).subs(symbols)
print "dx2[9][2][1] =",        cost.diff(vz3, vy2).subs(symbols)
print "dx2[9][2][2] =",        cost.diff(vz3, vz2).subs(symbols)
        
print "dx2[10][0][0] =",        cost.diff(vx3, vx3).subs(symbols)
print "dx2[10][0][1] =",        cost.diff(vx3, vy3).subs(symbols)
print "dx2[10][0][2] =",        cost.diff(vx3, vz3).subs(symbols)
print "dx2[10][1][0] =",        cost.diff(vy3, vx3).subs(symbols)
print "dx2[10][1][1] =",        cost.diff(vy3, vy3).subs(symbols)
print "dx2[10][1][2] =",        cost.diff(vy3, vz3).subs(symbols)
print "dx2[10][2][0] =",        cost.diff(vz3, vx3).subs(symbols)
print "dx2[10][2][1] =",        cost.diff(vz3, vy3).subs(symbols)
print "dx2[10][2][2] =",        cost.diff(vz3, vz3).subs(symbols)
        
print "dx2[11][0][0] =",        cost.diff(vx3, vx4).subs(symbols)
print "dx2[11][0][1] =",        cost.diff(vx3, vy4).subs(symbols)
print "dx2[11][0][2] =",        cost.diff(vx3, vz4).subs(symbols)
print "dx2[11][1][0] =",        cost.diff(vy3, vx4).subs(symbols)
print "dx2[11][1][1] =",        cost.diff(vy3, vy4).subs(symbols)
print "dx2[11][1][2] =",        cost.diff(vy3, vz4).subs(symbols)
print "dx2[11][2][0] =",        cost.diff(vz3, vx4).subs(symbols)
print "dx2[11][2][1] =",        cost.diff(vz3, vy4).subs(symbols)
print "dx2[11][2][2] =",        cost.diff(vz3, vz4).subs(symbols)
        
print "dx2[12][0][0] =",        cost.diff(vx4, vx1).subs(symbols)
print "dx2[12][0][1] =",        cost.diff(vx3, vy1).subs(symbols)
print "dx2[12][0][2] =",        cost.diff(vx3, vz1).subs(symbols)
print "dx2[12][1][0] =",        cost.diff(vy4, vx1).subs(symbols)
print "dx2[12][1][1] =",        cost.diff(vy3, vy1).subs(symbols)
print "dx2[12][1][2] =",        cost.diff(vy3, vz1).subs(symbols)
print "dx2[12][2][0] =",        cost.diff(vz4, vx1).subs(symbols)
print "dx2[12][2][1] =",        cost.diff(vz3, vy1).subs(symbols)
print "dx2[12][2][2] =",        cost.diff(vz3, vz1).subs(symbols)
        
print "dx2[13][0][0] =",        cost.diff(vx4, vx2).subs(symbols)
print "dx2[13][0][1] =",        cost.diff(vx4, vy2).subs(symbols)
print "dx2[13][0][2] =",        cost.diff(vx4, vz2).subs(symbols)
print "dx2[13][1][0] =",        cost.diff(vy4, vx2).subs(symbols)
print "dx2[13][1][1] =",        cost.diff(vy4, vy2).subs(symbols)
print "dx2[13][1][2] =",        cost.diff(vy4, vz2).subs(symbols)
print "dx2[13][2][0] =",        cost.diff(vz4, vx2).subs(symbols)
print "dx2[13][2][1] =",        cost.diff(vz4, vy2).subs(symbols)
print "dx2[13][2][2] =",        cost.diff(vz4, vz2).subs(symbols)
        
print "dx2[14][0][0] =",        cost.diff(vx4, vx3).subs(symbols)
print "dx2[14][0][1] =",        cost.diff(vx4, vy3).subs(symbols)
print "dx2[14][0][2] =",        cost.diff(vx4, vz3).subs(symbols)
print "dx2[14][1][0] =",        cost.diff(vy4, vx3).subs(symbols)
print "dx2[14][1][1] =",        cost.diff(vy4, vy3).subs(symbols)
print "dx2[14][1][2] =",        cost.diff(vy4, vz3).subs(symbols)
print "dx2[14][2][0] =",        cost.diff(vz4, vx3).subs(symbols)
print "dx2[14][2][1] =",        cost.diff(vz4, vy3).subs(symbols)
print "dx2[14][2][2] =",        cost.diff(vz4, vz3).subs(symbols)
        
print "dx2[15][0][0] =",        cost.diff(vx4, vx4).subs(symbols)
print "dx2[15][0][1] =",        cost.diff(vx4, vy4).subs(symbols)
print "dx2[15][0][2] =",        cost.diff(vx4, vz4).subs(symbols)
print "dx2[15][1][0] =",        cost.diff(vy4, vx4).subs(symbols)
print "dx2[15][1][1] =",        cost.diff(vy4, vy4).subs(symbols)
print "dx2[15][1][2] =",        cost.diff(vy4, vz4).subs(symbols)
print "dx2[15][2][0] =",        cost.diff(vz4, vx4).subs(symbols)
print "dx2[15][2][1] =",        cost.diff(vz4, vy4).subs(symbols)
print "dx2[15][2][2] =",        cost.diff(vz4, vz4).subs(symbols)
"""